######################
## Tables and Plots for Solar Siting
## CEC insolation
######################
## 'NoCO2' output is exclusive of C02 impacts (cannot be apportioned by county of sink)

require(tigris)
require(tidycensus)
require(sp)
require(rgdal)
require(sf)
require(raster)
require(stargazer)
require(readxl)
require(haven)
require(zoo) 
require(snow)
require(doSNOW)
require(parallel)
require(foreach)
require(iterators)
require(foreign)
require(abind)


## Select output with CO2 damages, or without (uncomment one) ##
    # CO2 = 'noCO2' 
    CO2 = 'withCO2' 

    
WORK = file.path('Z:/user/ajk41/Solar siting/Scripts/JAERE Code Repository/Data JAERE')
WORK.OUT = file.path(WORK, 'intermediate_output')

PLOT.OUT = file.path(WORK, 'final_output',CO2)
dir.create(PLOT.OUT, recursive=T)

using.crs = 3395  #--> mercator. flat. Ugh.


### Load dzip_tot (the zip-by-pollutant and zip averted damages) ###
dzip_poll = readRDS(file.path(WORK.OUT, 'Damages','CEC_to_damages_by_POLL.rds'))[,.(zip, POLL, total_damages)]
  if(CO2=='noCO2') dzip_poll = dzip_poll[POLL=='CO2',total_damages:=0]

  dzip_tot = as.data.table(dzip_poll %>% group_by(zip) %>% summarize(total_damages = sum(total_damages)) %>%
                             ungroup %>% mutate(POLL = 'total_damages'))
  dzip_tot = rbindlist(list(dzip_tot, dzip_poll[,colnames(dzip_tot),with=F]))[order(zip),]
  dzip_tot = dcast(dzip_tot, zip  ~ POLL, value.var='total_damages')

  
### Load "later years" coefficients
lateryears = readRDS(file.path(WORK.OUT, 'Damages', 'LaterYearsDamages.rds'))
  if('list' %in% class(lateryears)) lateryears = rbindlist(lateryears)
  if(CO2=='noCO2') lateryears$CO2 = 0
  lateryears[,total_damages:=list(NOx+SO2+PM25+CO2)]

  
  
### Generation by zip
fzip = fread(file.path(WORK, 'solargen_zip_hour_CEC.csv'),
             col.names = c('junk','zip',paste0('hr_',1:8760)))[,-1][,zip:=sprintf("%05d", zip)]

fzip[,total_generation_MWh:=apply(fzip[,2:NCOL(fzip)], MAR=1, sum)/1000]

fzip = merge(x = fzip, y=readRDS(file.path(WORK, 'TEZ_mapping_zip_to_NERCsub.rds')), by='zip', all.x=F, all.y=F)
WEST.INT = c('AZNM','CAMX','NWPP','RMPA'); TRE.INT = c('ERCT')
fzip[,INT:=ifelse(pcac%in%TRE.INT, 'TRE',
                  ifelse(pcac%in%WEST.INT, 'WEST', 'EAST'))]
setcolorder(fzip, c('zip','pcac','INT','total_generation_MWh'))




### Texas boot damages 
TBD = readRDS(file.path(WORK.OUT, 'Damages','TexasBootsFinal.rds')) 
  if(CO2=='noCO2') TBD[,CO2:=0]
TBD[,total_damages:=NOx+SO2+PM25+CO2]
TBD = TBD[,j = list(boot_total_damages = mean(total_damages, na.rm=T),
                    sd_total_damages = sd(total_damages, na.rm=T)), by = 'zip']



### WECC boot damages
WBD = rbindlist(lapply(list.files(file.path(WORK.OUT,'Damages'), pattern='^WECCBootDamages_[[:upper:]]{4}\\.rds$', full.names = T), function(x) rbindlist(readRDS(x)[[1]])))
if(CO2=='noCO2') WBD[,CO2:=0]
WBD[,total_damages:=NOx+SO2+PM25+CO2]
WBD = WBD[,j = list(boot_total_damages=mean(total_damages, na.rm=T),
                      sd_total_damages = sd(total_damages, na.rm=T)), by=c('zip')]

### Combine them (they are exclusive zips - none overlap)
BD = bind_rows(WBD, TBD) %>% dplyr::arrange(zip) %>% setDT()
rm(WBD);rm(TBD)






#### Put them all together merging on ZIP ####
#---> Note: Borenstein and Bushnell data is reset to 0 in JAERE upload. 

## Zip code map and join data:
zip = st_as_sf(readRDS(file.path(WORK, 'FIPS Shapefiles', "USA_Zip_Code_Boundaries_v104_zip_poly.rds"))) %>% 
        dplyr::select(zip = czip, NAME, STATE, POPULATION, POP_SQMI) %>% 
  left_join(., fzip[,.(zip, pcac, INT, total_generation_MWh)], by='zip') %>%
  left_join(., dzip_tot, by='zip') %>%
  left_join(., fread(file.path(WORK, 'zip_subsidies.csv'), stringsAsFactors=F) %>% dplyr::mutate(zip = sprintf("%05d", zip)) %>% dplyr::select(zip, subsidy_kwh =2), by=c('zip' = 'zip')) %>%
  left_join(., readRDS(file.path(WORK, 'zip_to_system_lambda_BorensteinBushnell.rds')), by=c('zip' = 'zip')) %>%
  left_join(., fread(file.path(WORK, 'NumberSystems.csv')) %>% dplyr::mutate(zip = sprintf("%05d", zipcode)) %>% dplyr::select(zip, num_systems), by=c('zip'='zip')) %>% replace_na(list(num_systems = 0)) %>%
  left_join(., BD, by=c('zip'='zip')) %>%
  left_join(lateryears, by=c('zip'='zip'), suffix = c('','_later')) %>%
  dplyr::mutate( NET  = total_damages - (1000*subsidy_kwh*total_generation_MWh),
                 NET_later  = total_damages_later - (1000*subsidy_kwh*total_generation_MWh),
                 total_damages_diffLater = total_damages_later - total_damages,
                 NET_diffLater = NET_later - NET) %>%
  left_join(readRDS(file.path(WORK.OUT, 'ShareCaptured.rds'))[,.(zip, pct_own_state)], by = c('zip' = 'zip') ) %>%
  st_transform(crs=using.crs) %>% 
  st_simplify(dTolerance = 500) %>% 
  st_buffer(dist=0)



# write out data for reallocation simulations #
write_csv(zip %>% st_set_geometry(NULL), file.path(PLOT.OUT, paste0('Generation_and_damages_',CO2,'_',Sys.Date(),'_eGrid_SR.csv')))




#### Table 1: Estimated Avoided Damages Summary Statistics ####
# Not in this script


#### Table 2:  ####
# Not in this script


#### Figure 1: Total Annual A/C Electricity Generation per 4kw Solar Capacity ####
# Figure 1: Annual AC Electricity Generation per 1kW Capacity Across the U.S.

pdf(file.path(PLOT.OUT, 'Solar_insolation_CEC_4kw.pdf'), width=8, height=4)
ggplot(zip, aes(col=total_generation_MWh, fill=total_generation_MWh), lwd=0) + geom_sf() +
  scale_fill_viridis_c(labels=prettyNum, breaks = seq(from=0, to=8, by=.5), option='D', guide=guide_colorbar(title='Total Annual MWh Generated\nper 4kW capacity', position='bottomleft')) +
  scale_color_viridis_c(breaks = seq(from=0, to=8, by=.5), option='D', guide=F) +
  theme(text = element_text(color = "black")
        ,plot.title = element_text(size = 18, color = "black")
        ,plot.background = element_blank()
        ,panel.background = element_blank()
        ,panel.grid = element_line(color='white')
        ,axis.text = element_blank()
        ,axis.ticks = element_blank()
        ,legend.background = element_blank()
  )
dev.off()





#### Figure 2: Damages Per Ton SO2 Emissions by County ####
AP3 = readRDS(file.path(WORK.OUT, 'AP3.rds')) # you need package abind loaded

APSO2 = apply(AP3[,,'SO2'], MAR=1, sum) %>%
  as.tibble() %>%
  dplyr::mutate(fips = dimnames(AP3)[[1]]) %>%
  dplyr::select(fips, APSO2_damages = 1) %>%
  dplyr::mutate(fips = replace(fips, fips=='12025', '12086')) # Miami added/renamed a county...

fipsco = st_read(dsn=file.path(WORK, 'FIPS Shapefiles'), layer='cb_2013_us_county_500k') %>%
  st_transform(crs = st_crs(zip)[[2]]) %>%
  dplyr::select(fips = GEOID, coname = NAME) %>%
  dplyr::mutate(fips = replace(fips, fips=='08014','08013')) %>% # Broomfield Co. formed from Boulder County (013) in 2001
  dplyr::mutate(fips = replace(fips, fips=='46102','46113')) %>% # Shannon Co. ND (46113) was renamed and renumbered Oglala Lakota (46102) in 2015
  left_join(APSO2, by=c('fips'='fips'))

write_csv(fipsco %>% dplyr::arrange(-APSO2_damages) %>% st_set_geometry(NULL),
          file.path(PLOT.OUT, 'Damages_per_1ton_SO2_in_fips.csv'))

## plot it - Figure SOX ##
pdf(file.path(PLOT.OUT, paste0('Damages_per_1ton_SO2_in_fips.pdf')), width=8, height=4)
ggplot(fipsco, aes(col=APSO2_damages, fill=APSO2_damages), lwd=0) + geom_sf() +
  scale_fill_viridis_c(labels=scales::dollar, option='D', guide=guide_colorbar(title='Damages per ton\nSO2 emitted', position='bottomleft')) + 
  scale_color_viridis_c(option='D', guide=F) +
  theme(text = element_text(color = "black")
        ,plot.title = element_text(size = 18, color = "black")
        ,plot.background = element_blank()
        ,panel.background = element_blank()
        ,panel.grid = element_line(color='white')
        ,axis.text = element_blank()
        ,axis.ticks = element_blank()
        ,legend.background = element_blank()
  )
dev.off()



#### Figure 3: Total Annual Avoided Environmental Damages per 4kw Solar Capacity ####
pdf(file.path(PLOT.OUT, paste0('Total_averted_damages_CEC_4kw_figure3_oldfig5_',CO2,'.pdf')), width=8, height=4)
ggplot(zip, aes(col=total_damages, fill=total_damages), lwd=0) + geom_sf() +
  scale_fill_viridis_c(labels = scales::dollar, option='D', guide=guide_colorbar(title='Total Damages Avoided\nper 4kW capacity', position='bottomleft')) + 
  scale_color_viridis_c(option='D', guide=F) +
  theme(text = element_text(color = "black")
        ,plot.title = element_text(size = 18, color = "black")
        ,plot.background = element_blank()
        ,panel.background = element_blank()
        ,panel.grid = element_line(color='white')
        ,axis.text = element_blank()
        ,axis.ticks = element_blank()
        ,legend.background = element_blank()
  )
dev.off()

#### Table 3: Simulated Allocations of Installed Solar Capacity ####
# Not in this script

#### Table 4: In-state vs Out-of-state Environmental Benefits ####
state_to_state_installed_cap_T4 = readRDS(file.path(WORK.OUT, 'Table4raw_state_to_state_installed_cap.rds')) # From state_to_state_installed_cap
# Formatting is left as an exercise for the reader

#### Table 5: Environmental Benefit Gain from out-of-state capacity relative to in-state capacity ####
state_to_state_installed_cap_T5 = readRDS(file.path(WORK.OUT, 'Table5raw_state_to_state_installed_cap.rds')) # From state_to_state_installed_cap
# Formatting is left as an exercise for the reader

#### Figure 4: Installed Rooftop Solar Arrays ####
# Not in this script

#### Figure 5: Intrastate Reallocation of Solar Capacity subject to Grid Stability Concerns ####
# Not in this script

#### Figure 6: Interstate Reallocation of Solar Capacity subject to Grid Stability Concerns ####
# Not in this script





##################
#### Appendix ####

#### Table A.1: Benefits to Largest Out-of-state Beneficiary ####
state_to_state_installed_capacity_A1 = readRDS(file.path(WORK.OUT, 'TableA1raw_state_to_state_installed_cap.rds')) # From state_to_state_installed_cap
# Formatting is left as an exercise for the reader


#### Table A.2: Energy Value and Total Benefits by State ####
A2 = zip %>% group_by(STATE) %>% st_set_geometry(NULL) %>% 
  dplyr::summarize(anntot = weighted.mean(total_damages, w = num_systems), annEV = weighted.mean(genval_severin, w = num_systems)) %>% 
  dplyr::mutate(tottot = round(anntot + annEV,0)) %>% 
  dplyr::mutate(anntot = round(anntot, 0), 
                annEV = round(annEV, 0)) %>%
  arrange(-tottot) %>%
  dplyr::select(State = STATE, AnnTotalBen = tottot, AnnEnergyVal = annEV)

stargazer(A2, type='latex',summary=F, digits = 0,
          out = file.path(PLOT.OUT, paste0('Table_A2_',CO2,'_',Sys.Date(),'.tex')))
stargazer(A2, type='html', summary=F, digits = 0,
          out = file.path(PLOT.OUT, paste0('Table_A2_',CO2,'_',Sys.Date(),'.html')))




#### Table A.3: State Average Historical Subsidies and Avoided Damages ####
# Not in this script

#### Table A.4: Simulated Allocation of Installed Solar Capacity to Maximize Environmental Plus Energy Benefits ####
# Not in this script

#### Table A.5: Uncertainty in Gains from Capacity Reallocation ####
# Not in this script

#### Figure A.1 & A.2: Marginal Emissions by Month and Daylight House and 95% CI ####
Blist_collapse = rbindlist(lapply(readRDS(WORK.OUT, 'BlistTexas.rds'), function(x){
  x %>% dplyr::filter(abs(GLOAD_MASS)<10 & !grepl('SOL_', PCA)) %>% # trim those with extremely bad fit (winsorized below)
    group_by(intx, PCA, mo, hour,iterB) %>%
    summarize_at(vars(ends_with('MASS')), sum, na.rm=T) %>%
    dplyr::arrange(mo, hour, iterB) %>% ungroup() %>%
    dplyr::mutate(hour = hour + 1) # 0-23 ---> 1-24 for he.
}))

use.vars = c('GLOAD_MASS','SO2_MASS','CO2_MASS','NOX_MASS')

for(plot_var in use.vars){
  print(plot_var)
  
  Blist_conf = Blist_collapse %>% dplyr::filter(hour>=6 & hour<=21) %>%
    group_by(mo, hour) %>% summarize(conf.low = quantile(get(plot_var), .025),
                                     conf.high = quantile(get(plot_var), .975))
  
  p = ggplot(Blist_collapse %>% dplyr::filter(hour>=6 & hour<=21), aes_string(x = 'hour', y = plot_var)) +
    geom_path(aes(group = iterB), col = alpha('gray', .2)) +
    geom_path(data = Blist_conf, aes(x = hour, y = conf.low), col='gray20') +
    geom_path(data = Blist_conf, aes(x = hour, y = conf.high), col='gray20') +
    geom_path(data = Blist_collapse %>% dplyr::filter(hour>=6 & hour<=21 & iterB==2), col='blue') + # the estimation code places the full (main) samle est. at iterB==2.
    scale_y_continuous(name = strsplit(plot_var, '_')[[1]][1]) +
    scale_x_continuous(name = 'Hour', breaks=seq(6, 21, by=4)) +
    facet_grid(cols=vars(mo), labeller = labeller(mo = function(x) {month.abb[as.integer(x)]})) + # must take char. vector, return char vector.
    theme_bw(base_size = 8)
  
  ggsave(p, filename = file.path(dirname(PLOT.OUT), paste0('A1A2',plot_var,'_plot.png')), width=6.5, height=3, dpi='print', units='in')
}

rm(Blist)
rm(Blist_conf)
rm(Blist_collapse)


#### Figure A.3: Annual Value of Electricity Generation per kW ####
# See Kirkpatrick (2018). Uses confidential pricing node locations. 

#### Table A.6: Sensitivity of Avoided Damages to Grid Evolution ####
stargazer(zip %>% st_set_geometry(NULL) %>% dplyr::select(total_damages_later, CO2_later, NOx_later, PM25_later, SO2_later), type='latex',summary=T, 
          out = file.path(PLOT.OUT, paste0('Table_1_',CO2,'_',Sys.Date(),'eGrid_SR.tex')))
stargazer(zip %>% st_set_geometry(NULL) %>% dplyr::select(total_damages_later, CO2_later, NOx_later, PM25_later, SO2_later), type='html',summary=T, 
          out = file.path(PLOT.OUT, paste0('Table_1_',CO2,'_',Sys.Date(),'eGrid_SR.html')))


#### Figure A.4: Interstate Reallocation of Solar Capacity for Environmental and Energy Benefits subject to Grid Stability Concerns  ####
# Not in this script

#### Figure A.5: Intrastate Reallocation of Solar Capacity for Environmental and Energy Benefits subject to Grid Stability Concerns ####
# Not in this script


#### Figure A.6: In-State Shares of Avoided Local Pollution Benefits ####
png(file.path(PLOT.OUT, paste0('benefits_in_out_state_noCO2.png')), width=8, height=4, res=300, units='in')
ggplot(zip, aes(col=pct_own_state, fill=pct_own_state), lwd=0) + geom_sf() +
  scale_fill_viridis_c(breaks = seq(from=0, to=1, by=.2), option='D', guide=guide_colorbar(title='Share of benefits\ncaptured in-state', position='bottomleft', labeller=scales::percent)) + 
  scale_color_viridis_c(option='D', guide=F) +
  theme(text = element_text(color = "black")
        ,plot.title = element_text(size = 18, color = "black")
        ,plot.background = element_blank()
        ,panel.background = element_blank()
        ,panel.grid = element_line(color='white')
        ,axis.text = element_blank()
        ,axis.ticks = element_blank()
        ,legend.background = element_blank()
  )
dev.off()

#### Figure A.7: Distribution of Reported System Azimuth in California ####
# Not in this script























#### Unused ####
#### Delete ####

# 
# 
# pdf(file.path(PLOT.OUT, paste0('Total_averted_damages_later_minus_total_damages_CEC_4kw_',CO2,'.pdf')), width=8, height=4)
# ggplot(zip, aes(col=total_damages_diffLater, fill=total_damages_diffLater), lwd=0) + geom_sf() +
#   scale_fill_viridis_c(labels = scales::dollar, option='D', guide=guide_colorbar(title='Difference in Damages Avoided\nper 4kW capacity', position='bottomleft')) + 
#   scale_color_viridis_c(option='D', guide=F) +
#   theme(text = element_text(color = "black")
#         ,plot.title = element_text(size = 18, color = "black")
#         ,plot.background = element_blank()
#         ,panel.background = element_blank()
#         ,panel.grid = element_line(color='white')
#         ,axis.text = element_blank()
#         ,axis.ticks = element_blank()
#         ,legend.background = element_blank()
#   )
# dev.off()
# 
# 
# 
# ## Not a figure, but do it anyways:
# pdf(file.path(PLOT.OUT, paste0('Total_averted_damages_CA_CEC_4kw_',CO2,'.pdf')), width=8, height=4)
# ggplot(zip %>% dplyr::filter(STATE=='CA'), aes(col=total_damages, fill=total_damages), lwd=0) + geom_sf() +
#   scale_fill_viridis_c(option='D', guide=guide_colorbar(title='Total Damages Avoided\nper 4kW capacity', position='bottomleft'), labels=scales::dollar) + 
#   scale_color_viridis_c(option='D', guide=F) +
#   theme(text = element_text(color = "black")
#         ,plot.title = element_text(size = 18, color = "black")
#         ,plot.background = element_blank()
#         ,panel.background = element_blank()
#         ,panel.grid = element_line(color='white')
#         ,axis.text = element_blank()
#         ,axis.ticks = element_blank()
#         ,legend.background = element_blank()
#   )
# dev.off()
# 
# 
# ### Figure 6:
# pdf(file.path(PLOT.OUT, paste0('CO2_averted_damages_CEC_4kw_figure6_',CO2,'.pdf')), width=8, height=4)
# ggplot(zip, aes(col=CO2, fill=CO2), lwd=0) + geom_sf() +
#   scale_fill_viridis_c(option='D', guide=guide_colorbar(title='Total CO2 Damages Avoided\nper 4kW capacity', position='bottomleft'), labels=scales::dollar) + 
#   scale_color_viridis_c(option='D', guide=F) +
#   theme(text = element_text(color = "black")
#         ,plot.title = element_text(size = 18, color = "black")
#         ,plot.background = element_blank()
#         ,panel.background = element_blank()
#         ,panel.grid = element_line(color='white')
#         ,axis.text = element_blank()
#         ,axis.ticks = element_blank()
#         ,legend.background = element_blank()
#   )
# dev.off()
# 
# 
# ### Figure 7:
# pdf(file.path(PLOT.OUT, 'SO2_averted_damages_CEC_4kw_figure7.pdf'), width=8, height=4)
# ggplot(zip, aes(col=SO2, fill=SO2), lwd=0) + geom_sf() +
#   scale_fill_viridis_c(breaks = seq(from=0, to=4200, by=200), option='D', guide=guide_colorbar(title='SO2 Damages Avoided\nper 4kW capacity', position='bottomleft'), labels=scales::dollar) + 
#   scale_color_viridis_c(option='D', guide=F) +
#   theme(text = element_text(color = "black")
#         ,plot.title = element_text(size = 18, color = "black")
#         ,plot.background = element_blank()
#         ,panel.background = element_blank()
#         ,panel.grid = element_line(color='white')
#         ,axis.text = element_blank()
#         ,axis.ticks = element_blank()
#         ,legend.background = element_blank()
#   )
# dev.off()
# 
# 
# ### Figure 8:
# pdf(file.path(PLOT.OUT, 'NOX_averted_damages_CEC_4kw_figure8.pdf'), width=8, height=4)
# ggplot(zip, aes(col=NOx, fill=NOx), lwd=0) + geom_sf() +
#   scale_fill_viridis_c(breaks = seq(from=0, to=1700, by=20), option='D', guide=guide_colorbar(title='NOX Damages Avoided\nper 4kW capacity', position='bottomleft'), labels=scales::dollar) + 
#   scale_color_viridis_c(option='D', guide=F) +
#   theme(text = element_text(color = "black")
#         ,plot.title = element_text(size = 18, color = "black")
#         ,plot.background = element_blank()
#         ,panel.background = element_blank()
#         ,panel.grid = element_line(color='white')
#         ,axis.text = element_blank()
#         ,axis.ticks = element_blank()
#         ,legend.background = element_blank()
#   )
# dev.off()
# 
# ### Figure 9:
# pdf(file.path(PLOT.OUT, 'PM25_averted_damages_CEC_4kw_figure9.pdf'), width=8, height=4)
# ggplot(zip, aes(col=PM25, fill=PM25), lwd=0) + geom_sf() +
#   scale_fill_viridis_c(breaks = seq(from=0, to=250, by=10), option='D', guide=guide_colorbar(title='PM2.5 Damages Avoided\nper 4kW capacity', position='bottomleft'), labels=scales::dollar) + 
#   scale_color_viridis_c(option='D', guide=F) +
#   theme(text = element_text(color = "black")
#         ,plot.title = element_text(size = 18, color = "black")
#         ,plot.background = element_blank()
#         ,panel.background = element_blank()
#         ,panel.grid = element_line(color='white')
#         ,axis.text = element_blank()
#         ,axis.ticks = element_blank()
#         ,legend.background = element_blank()
#   )
# dev.off()
# 
# 
# ### Figure 10:
# pdf(file.path(PLOT.OUT, paste0('Net_averted_damages_CEC_4kw_figure10_',CO2,'.pdf')), width=8, height=4)
# ggplot(zip, aes(col=NET, fill=NET), lwd=0) + geom_sf() +
#   scale_fill_viridis_c(breaks = c(-2000,-1500,-1000,-500,0,500), labels = scales::dollar,  option='D', guide=guide_colorbar(title='Damages Avoided Net of Incentives\nper 4kW capacity', position='bottomleft')) + 
#   scale_color_viridis_c(option='D', guide=F) +
#   theme(text = element_text(color = "black")
#         ,plot.title = element_text(size = 18, color = "black")
#         ,plot.background = element_blank()
#         ,panel.background = element_blank()
#         ,panel.grid = element_line(color='white')
#         ,axis.text = element_blank()
#         ,axis.ticks = element_blank()
#         ,legend.background = element_blank()
#   )
# dev.off()
# 
# 
# # 
# # ### Extra: pcac (eGrid subregions) by zip ###
# zip2 = zip %>%  st_set_precision(1000000) %>% lwgeom::st_make_valid() %>% group_by(pcac) %>% dplyr::summarize(.) %>% st_buffer(dist=1000)
# pdf(file.path(PLOT.OUT, 'combined_eGrid_subregions.pdf'), width=8, height=4)
# ggplot(zip2, aes(fill=pcac, col=pcac), lwd=0) + geom_sf() +
#   scale_fill_viridis_d(labels=prettyNum, option='D', guide=F) +
#   scale_color_viridis_d(option='D', guide=F) +
#   theme(text = element_text(color = "black")
#        ,plot.title = element_text(size = 18, color = "black")
#        ,plot.background = element_blank()
#        ,panel.background = element_blank()
#        ,panel.grid = element_line(color='white')
#        ,axis.text = element_blank()
#        ,axis.ticks = element_blank()
#        ,legend.background = element_blank())
# dev.off()
# 
# 
# 
# ## Figure ?? - System Lambda values - CO2 doesn't differ ##
# if(CO2=='withCO2'){
# pdf(file.path(PLOT.OUT, 'system_lambdas_severin.pdf'), width=8, height=4)
# ggplot(zip, aes(col=genval_severin, fill=genval_severin), lwd=0) + geom_sf() +
#   scale_fill_viridis_c(labels=scales::dollar, option='D', guide=guide_colorbar(title='System Lambda Lalue\nper 4kW capacity', position='bottomleft')) + 
#   scale_color_viridis_c(option='D', guide=F) +
#   theme(text = element_text(color = "black")
#         ,plot.title = element_text(size = 18, color = "black")
#         ,plot.background = element_blank()
#         ,panel.background = element_blank()
#         ,panel.grid = element_line(color='white')
#         ,axis.text = element_blank()
#         ,axis.ticks = element_blank()
#         ,legend.background = element_blank()
#   )
# dev.off()
# } # withCO2 only
# 
# 
# ## Texas (and WECC) SD of estimates ##
# if(CO2=='withCO2'){
# pdf(file.path(PLOT.OUT, 'TexasWECC_Bootstrap_withCO2only.pdf'), width=8, height=4)
# ggplot(zip %>% dplyr::filter(!is.na(sd_total_damages)), aes(col=sd_total_damages, fill=sd_total_damages), lwd=0) + geom_sf() +
#   scale_fill_viridis_c(labels=scales::dollar, option='D', guide=guide_colorbar(title='Std. Dev of Damages\nper 4kW capacity', position='bottomleft')) + 
#   scale_color_viridis_c(option='D', guide=F) +
#   theme(text = element_text(color = "black")
#         ,plot.title = element_text(size = 18, color = "black")
#         ,plot.background = element_blank()
#         ,panel.background = element_blank()
#         ,panel.grid = element_line(color='white')
#         ,axis.text = element_blank()
#         ,axis.ticks = element_blank()
#         ,legend.background = element_blank()
#   )
# dev.off()
# } # withCO2 only!
# 



# 
# ### Applying and plotting LaterYears ####
# ### Redo Table 1 with later results ###
# stargazer(zip %>% st_set_geometry(NULL) %>% dplyr::select(total_damages_later, NET_later, CO2_later, NOx_later, PM25_later, SO2_later), type='latex',summary=T, 
#           out = file.path(PLOT.OUT, paste0('Table_1_later_',CO2,'_',Sys.Date(),'eGrid_SR.tex')))
# stargazer(zip %>% st_set_geometry(NULL) %>% dplyr::select(total_damages_later, NET_later, CO2_later, NOx_later, PM25_later, SO2_later), type='html',summary=T, 
#           out = file.path(PLOT.OUT, paste0('Table_1_later_',CO2,'_',Sys.Date(),'eGrid_SR.html')))
# 
# ### Figure 5 Later:
# pdf(file.path(PLOT.OUT, paste0('Total_averted_damages_CEC_4kw_figure5_later_',CO2,'.pdf')), width=8, height=4)
# ggplot(zip, aes(col=total_damages_later, fill=total_damages_later), lwd=0) + geom_sf() +
#   scale_fill_viridis_c(labels = scales::dollar, option='D', guide=guide_colorbar(title='Total Damages Avoided\nper 4kW capacity', position='bottomleft')) + 
#   scale_color_viridis_c(option='D', guide=F) +
#   theme(text = element_text(color = "black")
#         ,plot.title = element_text(size = 18, color = "black")
#         ,plot.background = element_blank()
#         ,panel.background = element_blank()
#         ,panel.grid = element_line(color='white')
#         ,axis.text = element_blank()
#         ,axis.ticks = element_blank()
#         ,legend.background = element_blank()
#   )
# dev.off()
# 
# 
# ## Not a figure, but do it anyways_later:
# pdf(file.path(PLOT.OUT, paste0('Total_averted_damages_CA_CEC_4kw__later',CO2,'.pdf')), width=8, height=4)
# ggplot(zip %>% dplyr::filter(STATE=='CA'), aes(col=total_damages_later, fill=total_damages_later), lwd=0) + geom_sf() +
#   scale_fill_viridis_c(option='D', guide=guide_colorbar(title='Total Damages Avoided\nper 4kW capacity', position='bottomleft'), labels=scales::dollar) + 
#   scale_color_viridis_c(option='D', guide=F) +
#   theme(text = element_text(color = "black")
#         ,plot.title = element_text(size = 18, color = "black")
#         ,plot.background = element_blank()
#         ,panel.background = element_blank()
#         ,panel.grid = element_line(color='white')
#         ,axis.text = element_blank()
#         ,axis.ticks = element_blank()
#         ,legend.background = element_blank()
#   )
# dev.off()
# 
# 
# ### Figure 6:
# pdf(file.path(PLOT.OUT, paste0('CO2_averted_damages_CEC_4kw_figure6__later',CO2,'.pdf')), width=8, height=4)
# ggplot(zip, aes(col=CO2_later, fill=CO2_later), lwd=0) + geom_sf() +
#   scale_fill_viridis_c(option='D', guide=guide_colorbar(title='Total CO2 Damages Avoided\nper 4kW capacity', position='bottomleft'), labels=scales::dollar) + 
#   scale_color_viridis_c(option='D', guide=F) +
#   theme(text = element_text(color = "black")
#         ,plot.title = element_text(size = 18, color = "black")
#         ,plot.background = element_blank()
#         ,panel.background = element_blank()
#         ,panel.grid = element_line(color='white')
#         ,axis.text = element_blank()
#         ,axis.ticks = element_blank()
#         ,legend.background = element_blank()
#   )
# dev.off()
# 
# 
# ### Figure 7_later:
# pdf(file.path(PLOT.OUT, 'SO2_averted_damages_CEC_4kw_figure7_later.pdf'), width=8, height=4)
# ggplot(zip, aes(col=SO2_later, fill=SO2_later), lwd=0) + geom_sf() +
#   scale_fill_viridis_c(breaks = seq(from=0, to=4200, by=200), option='D', guide=guide_colorbar(title='SO2 Damages Avoided\nper 4kW capacity', position='bottomleft'), labels=scales::dollar) + 
#   scale_color_viridis_c(option='D', guide=F) +
#   theme(text = element_text(color = "black")
#         ,plot.title = element_text(size = 18, color = "black")
#         ,plot.background = element_blank()
#         ,panel.background = element_blank()
#         ,panel.grid = element_line(color='white')
#         ,axis.text = element_blank()
#         ,axis.ticks = element_blank()
#         ,legend.background = element_blank()
#   )
# dev.off()
# 
# 
# ### Figure 8_later:
# pdf(file.path(PLOT.OUT, 'NOX_averted_damages_CEC_4kw_figure8_later.pdf'), width=8, height=4)
# ggplot(zip, aes(col=NOx_later, fill=NOx_later), lwd=0) + geom_sf() +
#   scale_fill_viridis_c(breaks = seq(from=0, to=1700, by=20), option='D', guide=guide_colorbar(title='NOX Damages Avoided\nper 4kW capacity', position='bottomleft'), labels=scales::dollar) + 
#   scale_color_viridis_c(option='D', guide=F) +
#   theme(text = element_text(color = "black")
#         ,plot.title = element_text(size = 18, color = "black")
#         ,plot.background = element_blank()
#         ,panel.background = element_blank()
#         ,panel.grid = element_line(color='white')
#         ,axis.text = element_blank()
#         ,axis.ticks = element_blank()
#         ,legend.background = element_blank()
#   )
# dev.off()
# 
# ### Figure 9_later:
# pdf(file.path(PLOT.OUT, 'PM25_averted_damages_CEC_4kw_figure9_later.pdf'), width=8, height=4)
# ggplot(zip, aes(col=PM25_later, fill=PM25_later), lwd=0) + geom_sf() +
#   scale_fill_viridis_c(breaks = seq(from=0, to=250, by=10), option='D', guide=guide_colorbar(title='PM2.5 Damages Avoided\nper 4kW capacity', position='bottomleft'), labels=scales::dollar) + 
#   scale_color_viridis_c(option='D', guide=F) +
#   theme(text = element_text(color = "black")
#         ,plot.title = element_text(size = 18, color = "black")
#         ,plot.background = element_blank()
#         ,panel.background = element_blank()
#         ,panel.grid = element_line(color='white')
#         ,axis.text = element_blank()
#         ,axis.ticks = element_blank()
#         ,legend.background = element_blank()
#   )
# dev.off()
# 
# 
# ### Figure 10_later:
# pdf(file.path(PLOT.OUT, paste0('Net_averted_damages_CEC_4kw_figure10_later_',CO2,'.pdf')), width=8, height=4)
# ggplot(zip, aes(col=NET_later, fill=NET_later), lwd=0) + geom_sf() +
#   scale_fill_viridis_c(breaks = c(-2000,-1500,-1000,-500,0,500), labels = scales::dollar,  option='D', guide=guide_colorbar(title='Damages Avoided Net of Incentives\nper 4kW capacity', position='bottomleft')) + 
#   scale_color_viridis_c(option='D', guide=F) +
#   theme(text = element_text(color = "black")
#         ,plot.title = element_text(size = 18, color = "black")
#         ,plot.background = element_blank()
#         ,panel.background = element_blank()
#         ,panel.grid = element_line(color='white')
#         ,axis.text = element_blank()
#         ,axis.ticks = element_blank()
#         ,legend.background = element_blank()
#   )
# dev.off()
# 
# 
# ### 12-2019 draft Figure 7: Gain in local envrionmental benefits from optimal out-of-state capacity...
# pdf(file.path(PLOT.OUT, paste0('benefits_in_out_state_',CO2,'.pdf')), width=8, height=4)
# ggplot(zip, aes(col=pct_benefits_gain/100, fill=pct_benefits_gain/100), lwd=0) + geom_sf() +
#   scale_fill_viridis_c(labels = scales::percent, option='D', guide=guide_colorbar(title='Gain in appropriated\nbenefits per\n4kW capacity', position='bottomleft')) + 
#   scale_color_viridis_c(option='D', guide=F) +
#   theme(text = element_text(color = "black")
#         ,plot.title = element_text(size = 18, color = "black")
#         ,plot.background = element_blank()
#         ,panel.background = element_blank()
#         ,panel.grid = element_line(color='white')
#         ,axis.text = element_blank()
#         ,axis.ticks = element_blank()
#         ,legend.background = element_blank()
#   )
# dev.off()
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# ###############
# ### Here's some other processing stuff not used in the output anymore:
# 
# ## Texas data output for line 283 of paper, on the bootstrap results  ##
# # Replaced by WECC+Texas
# Blist = readRDS(jLoad(WORK.OUT, 'Blist'))
# Blist = Blist[which(!unlist(lapply(Blist, is.logical)))]  # only 752 got done out of 2000, the other 1250 are just logical NA's.
# 
# xx = Blist %>% bind_rows() %>%
#   dplyr::filter(!grepl('SOL_',PCA) & isSOL == 'qr.SOL' & hour>=6 & hour<=21) %>%
#   group_by(iterB, mo, hour) %>%
#   dplyr::summarize(GLOAD_MASS = sum(GLOAD_MASS, na.rm=T)) %>% 
#   ungroup() %>% group_by(iterB) %>%
#   dplyr::summarize(low.95 = quantile(GLOAD_MASS, .025),
#                    high.95= quantile(GLOAD_MASS, .975),
#                    mean.G = mean(GLOAD_MASS)) %>%
#   ungroup() %>%
#   dplyr::summarize_all(mean)
# 
# Sensitivity = Blist %>%
#   bind_rows() %>%
#   dplyr::filter(!grepl('SOL_',PCA) & isSOL == 'qr.SOL' & hour>=6 & hour<=21) %>%
#   group_by(iterB, mo, hour) %>%
#   dplyr::summarize_at(vars(contains('MASS')), sum) %>% 
#   ungroup() %>%
#   group_by(mo, hour) %>%
#   dplyr::summarize_at(.vars = vars(contains('MASS')), function(x) {
#     sd(x)/mean(x)}) %>%
#   ungroup() %>%
#   dplyr::summarize_at(.vars = vars(contains('MASS')), mean)
# Sensitivity
# 
# 
# 
# plant.coef = readRDS(file.path(BASE, 'ajk41/Solar Siting/Data/OUT_2018-10-03/betas_processed_v14.rds'))
# use.ORISPL = unique(plant.coef$ORISPL)   #  use.ORISPL = unique(as.numeric(sapply(strsplit(plant.coef$ORISPL, '_'), '[[', 2)))
# plant.coef[,ORISPL_BU:=ORISPL]
# #--? gets PM2.5 factor from old data.
# plant.coef = merge(x=plant.coef[,ORISPL:=paste0('ORISPL_', ORISPL)], 
#                    y = !!XX!! PM25,
#                    by = 'ORISPL')[,PM2_MASS:=GLOAD_MASS*PM2_factor][, c('ORISPL','PM2_factor') :=list(ORISPL, NULL)] 
# plant.coef.SOL   = plant.coef %>% dplyr::filter(!grepl('SOL_', PCA) & isSOL=='qr.SOL')  %>% setDT()
# 
# 
# 
# 

